Cylindrical Kink

M. E. Mauel
Dept. of Applied Physics
Columbia University

Introduction

Off[Power::infy]

This notebook uses NDSolve[...] to compute the external kink eigenfunction for the Wesson current profile.

Alfvén Frequency

The Alfvén frequency is given by ((Subscript[ω, A])^2) = ((((Subscript[B, t])/(Subscript[q, a])R))^2)/(Subscript[μ, 0])ρ. 

ωA1 is the frequency that corresponds to (Subscript[q, a]) = 1.

Plasma Resistivity (MKS units : Ohm meter) ...

	η = 1.03*^-4 1.5 14 / ((Subscript[T, e])^(3/2)) with (Subscript[T, e]) in eV

ηHBTEP = 1.03*^-4 1.5 14/(temp^(3/2)) /. temp -> 25.0;
% // EngineeringForm

hbtepRule = {b0 -> 0.35, ρ -> N[10.0*^12 1.0*^6 * 2 * 1.67*^-27], R0 -> 2 π, a0 -> 0.14, μ0 -> 4 π (10^(-7)),η -> ηHBTEP}

ωA1 = Sqrt[(b0/R0)^2/μ0/ρ ]/. hbtepRule (* ωA((Subscript[q, a])= 3) *);
ωA1/2/π // EngineeringForm

Current Density (Wesson)

jz0[x_] := 2 (b0/μ0/R0/q0) (1 - (x^2))^νj

(a0^2) Integrate[2 π x jz0[x], {x, 0, 1}, Assumptions -> {νj > 1}]
2 π (a0^2) b0/(μ0 R0 %) // Simplify (* q(a) *)

Safety Factor

q[x_] := qa (x^2)/ (1 - (1 - (x^2))^(νj + 1))

Plot[{q[x] /. qa -> 2.9 /. νj -> 1.1,q[x] /. qa -> 2.9 /. νj -> 2,q[x] /. qa -> 2.9 /. νj -> 4}, {x, 0.001, 1}, Frame -> True, FrameLabel -> {"r/a","q(r)"},PlotRange -> {0,3},PlotLabel -> "Wesson Safety Factor for νj = 1.1, 2.0, 4.0", GridLines -> {None, {1,2}}]

Plot[{jz0[x](10^(-6)) /. q0 -> qa/(1 + νj)/. qa -> 2.9 /. νj -> 1.1 /. hbtepRule,jz0[x](10^(-6)) /. q0 -> qa/(1 + νj)/. qa -> 2.9 /. νj -> 2.0 /. hbtepRule,jz0[x](10^(-6)) /. q0 -> qa/(1 + νj)/. qa -> 2.9 /. νj -> 4.0 /. hbtepRule}, {x, 0, 1}, Frame -> True, FrameLabel -> {"r/a","Jz(r) (MA/(m^2))"},PlotLabel -> "Wesson Current Density for νj = 1.1, 2.0, 4.0"]

Vacuum Boundary Condition

The perturbed helical flux between the plasma (r = a) and an ideal wall (r = b) is

ψVac[r_] := ψa ((b/r)^m - (r/b)^m)/((b/a)^m - (a/b)^m)

Plot[{ψVac[r]/ψa /. b-> 1.8 /. a -> 1 /. m -> 3,ψVac[r]/ψa /. b-> 1.8 /. a -> 1 /. m -> 5}, {r, 1, 1.8}, PlotLabel -> "Vacuum Perturbed Flux"]

D[ψVac[r],r]/ψa // Simplify

% /. r -> a // Simplify

which we can write as

- (m/a) ( (1 + (a/b)^(2m)) / (1 - (a/b)^(2m))
or
- (m/a) (1 + c)/(1-c), where c = (a/b)^(2m) is the “coupling coefficient” between the plasma boundary and the wall.

ψpψBoundary[a_,b_,m_] := - (m/a) ((1 + (((a/b))^(2m)))/(1 - (((a/b))^(2m))))

Eqs for MHD Force Balance

With x = r/a and the Wesson current profile, the equation for MHD force balance is

eqsIdeal = {D[ψ[x],x] == ψp[x], x D[x ψp[x], x] - m^2 ψ[x] + 4 νj (m/n/q0) x^2 (1 - x^2)^(νj -1) ψ[x]/(m/n/q[x] - 1) == 0}

Example Solutions

Because of (integrable/resolvable) singularities in q[x] at the origin and at the edge, the eigenfunction is integrated up-to but not including these limits.

δ0 = 0.0001; (* how close to integrate to boundaries *)

origin = {ψ[δ0] == 0.0, ψp[δ0] == 1.0}

((Subscript[q, a]), (Subscript[q, 0])) = (2.8, 1.8)

gqa = 2.8;
gq0 = 1.8;

sol = NDSolve[Evaluate[eqsIdeal ~Join~ origin /. νj -> qa/q0 - 1 /. qa -> gqa/. q0 -> gq0 /. m -> 3 /. n -> 1], {ψ[x], ψp[x]}, {x, δ0, 1.0-δ0}] // Flatten

pp=Plot[Evaluate[((ψ[x]/.sol/.x->xx)/(ψ[x]/.sol/.x->1-δ0))],{xx,δ0,1-δ0},PlotLabel->"ψ[x]"]

Vacuum Region (Ideal Wall)

Λ[b_,m_] := (1 + (1/b)^(2m)) / (1 - (1/b)^(2m))

ψVac[b_,m_][x_] := (* normalized to 1 *) ( (b/x)^m - (x/b)^m )/( b^m - (1/b)^m )

pv=Plot[{ψVac[1.8,3][x],If[x < 1.1,ψVac[1.1,3][x],0.0]},{x,1,1.8}, PlotLabel -> "Perturbed Vacuum Flux with b/a = 1.8", PlotRange -> {0,1}]

Show[pp,pv,PlotLabel->"ψ[x] ((Subscript[q, a]), (Subscript[q, 0])) = "<>ToString[{gqa,gq0}], PlotRange -> {{0,1.5}, Automatic}, Frame -> True, FrameLabel -> {"r/a", "ψ"}]

Growth Rate

Δa = (D[ψVac[1.8,3][x],x] /. x -> 1) - (ψp[x]/ψ[x] /. sol /. x -> 1 - δ0)

(ψp[x]/ψ[x] /. sol /. x -> 1-δ0) + 2 m/(m - n qa) /. m -> 3 /. n -> 1 /. qa -> gqa

(* Frequency ω (unstable with I) *) Sqrt[(-(ωA1/qa)^2 (((m - n qa))^2) Δa (1/%) /. qa -> gqa/. q0 -> gq0 /. m -> 3 /. n -> 1)]

Stream Function

χs[x_] = - ψ[x]/(m/q[x] - n)

χsp[x_] = D[χs[x],x] /. D[ψ[x], x] -> ψp[x]

χsEdge = χs[x] /. sol  /. νj -> qa/q0 - 1/. qa -> gqa /. q0 -> gq0/. m -> 3 /. n -> 1 /. x -> 1 - δ0

Plot[Evaluate[((χs[x])/(χsEdge))/.sol/.νj->(qa/q0)-1/.qa->gqa/.q0->gq0/.m->3/.n->1],{x,δ0,1-δ0},PlotRange->All,PlotLabel->"(χ[x]/ω) ((Subscript[B, t])/(Subscript[R, 0])) Normalized"]

Energy Integrals

Kinetic Term, δK...

(1/2)ρ ( D[χs[x],x]^2 + (χs[x]/x)^2) (b0/R0)^(-2)

(1/2)ρ ( χsp[x]^2 + (χs[x]/x)^2) (b0/R0)^(-2) /. sol  /. νj -> qa/q0 - 1/. qa -> gqa /. q0 -> gq0/. m -> 3 /. n -> 1

(0.5 ρ )(b0/R0)^(-2) (qa/3)^2 NIntegrate[Evaluate[( χsp[x]^2 + (χs[x]/x)^2) /. sol/. νj -> qa/q0 - 1/. qa -> gqa /. q0 -> gq0/. m -> 3 /. n -> 1], {x, δ0, 1 - δ0}]/(π m ψ[x]^2/μ0  /. sol/. m -> 3 /. x -> 1 - δ0)

((Subscript[q, a]), (Subscript[q, 0])) = (2.8, 1.1)

gqa = 2.8;
gq0 = 1.1;

sol = NDSolve[Evaluate[eqsIdeal ~Join~ origin /. νj -> qa/q0 - 1 /. qa -> gqa/. q0 -> gq0 /. m -> 3 /. n -> 1], {ψ[x], ψp[x]}, {x, δ0, 1.0-δ0}] // Flatten

pp=Plot[Evaluate[((ψ[x]/.sol/.x->xx)/(ψ[x]/.sol/.x->1-δ0))],{xx,δ0,1-δ0},PlotLabel->"ψ[x]"]

Vacuum Region (Ideal Wall)

Λ[b_,m_] := (1 + (1/b)^(2m)) / (1 - (1/b)^(2m))

ψVac[b_,m_][x_] := (* normalized to 1 *) ( (b/x)^m - (x/b)^m )/( b^m - (1/b)^m )

pv=Plot[ψVac[1.8,3][x],{x,1,1.8}]

Show[pp,pv,PlotLabel->"ψ[x] ((Subscript[q, a]), (Subscript[q, 0])) = "<>ToString[{gqa,gq0}], PlotRange -> {{0,1.5}, Automatic}, Frame -> True, FrameLabel -> {"r/a", "ψ"}]

Growth Rate

Δa = (D[ψVac[1.8,3][x],x] /. x -> 1) - (ψp[x]/ψ[x] /. sol /. x -> 1 - δ0)

(ψp[x]/ψ[x] /. sol /. x -> 1-δ0) + 2 m/(m - n qa) /. m -> 3 /. n -> 1 /. qa -> gqa

(* Frequency ω (unstable with I) *) Sqrt[(-(ωA1/qa)^2 (((m - n qa))^2) Δa (1/%) /. qa -> gqa/. q0 -> gq0 /. m -> 3 /. n -> 1)]

((Subscript[q, a]), (Subscript[q, 0])) = (2.8, 0.9)

gqa = 2.8;
gq0 = 0.9;

sol = NDSolve[Evaluate[eqsIdeal ~Join~ origin /. νj -> qa/q0 - 1 /. qa -> gqa/. q0 -> gq0 /. m -> 3 /. n -> 1], {ψ[x], ψp[x]}, {x, δ0, 1.0-δ0}] // Flatten

pp=Plot[Evaluate[((ψ[x]/.sol/.x->xx)/(ψ[x]/.sol/.x->1-δ0))],{xx,δ0,1-δ0},PlotLabel->"ψ[x]"]

Vacuum Region (Ideal Wall)

Λ[b_,m_] := (1 + (1/b)^(2m)) / (1 - (1/b)^(2m))

ψVac[b_,m_][x_] := (* normalized to 1 *) ( (b/x)^m - (x/b)^m )/( b^m - (1/b)^m )

pv=Plot[ψVac[1.8`,3][x],{x,1,1.8`}]

Show[pp,pv,PlotLabel->"ψ[x] ((Subscript[q, a]), (Subscript[q, 0])) = "<>ToString[{gqa,gq0}], PlotRange -> {{0,1.5}, Automatic}, Frame -> True, FrameLabel -> {"r/a", "ψ"}]

Growth Rate

Δa = (D[ψVac[1.8,3][x],x] /. x -> 1) - (ψp[x]/ψ[x] /. sol /. x -> 1 - δ0)

(ψp[x]/ψ[x] /. sol /. x -> 1-δ0) + 2 m/(m - n qa) /. m -> 3 /. n -> 1 /. qa -> gqa

(* Frequency ω (unstable with I) *) Sqrt[(-(ωA1/qa)^2 (((m - n qa))^2) Δa (1/%) /. qa -> gqa/. q0 -> gq0 /. m -> 3 /. n -> 1)]

Eigenvalue Function

ωMHD[qqa_,qq0_,bwall_:1.8,mm_:3,nn_:1]:= Module[{sol,Δa,Δχa},
sol = NDSolve[Evaluate[eqsIdeal ~Join~ origin /. νj -> qa/q0 - 1 /. qa -> qqa/. q0 -> qq0 /. m -> mm /. n -> nn], {ψ[x], ψp[x]}, {x, δ0, 1.0-δ0}] // Flatten;
Δa = (D[ψVac[bwall,3][x],x] /. x -> 1) - (ψp[x]/ψ[x] /. sol /. x -> 1 - δ0);
Δχa = (ψp[x]/ψ[x] /. sol /. x -> 1-δ0) + 2 m/(m - n qa) /. m -> mm /. n -> nn /. qa -> qqa;
Sqrt[(-(ωA1/qa)^2 (((m - n qa))^2) Δa (1/Δχa) /. qa -> qqa/. q0 -> qq0 /. m -> mm /. n -> nn)]]

ωMHD[2.9,1.1]
ωMHD[2.9,1.8]

Plot[{Re[ωMHD[qa,1.5]],Im[ωMHD[qa,1.5]]}, {qa, 2.5, 3.2}];
Plot[{Re[ωMHD[qa,1.4]],Im[ωMHD[qa,1.4]]}, {qa, 2.5, 3.2}];
Plot[{Re[ωMHD[qa,1.3]],Im[ωMHD[qa,1.3]]}, {qa, 2.5, 3.2}];
Show[%%%,%%,%, Frame -> True, FrameLabel -> {"Edge Safety Factor ((Subscript[q, a]))","ω, γ"}, PlotLabel -> "External Kink for Wesson Profile, (Subscript[q, 0]) = (1.5, 1.4, 1.3)", PlotRange -> {0,20000}]

Plot[{Re[ωMHD[qa,1.5,1.2]],Im[ωMHD[qa,1.5,1.2]]}, {qa, 2.5, 3.2}];
Plot[{Re[ωMHD[qa,1.4,1.2]],Im[ωMHD[qa,1.4,1.2]]}, {qa, 2.5, 3.2}];
Plot[{Re[ωMHD[qa,1.3,1.2]],Im[ωMHD[qa,1.3,1.2]]}, {qa, 2.5, 3.2}];
Show[%%%,%%,%, Frame -> True, FrameLabel -> {"Edge Safety Factor ((Subscript[q, a]))","ω, γ"}, PlotLabel -> "External Kink for Wesson Profile, (Subscript[q, 0]) = (1.5, 1.4, 1.3)", PlotRange -> {0,20000}]

Plot[{Re[ωMHD[qa,1.5,1.8,5,2]],Im[ωMHD[qa,1.5,1.8,5,2]]}, {qa, 2.1, 3.2}];
Plot[{Re[ωMHD[qa,1.4,1.8,5,2]],Im[ωMHD[qa,1.4,1.8,5,2]]}, {qa, 2.1, 3.2}];
Plot[{Re[ωMHD[qa,1.3,1.8,5,2]],Im[ωMHD[qa,1.3,1.8,5,2]]}, {qa, 2.1, 3.2}];
Show[%%%,%%,%, Frame -> True, FrameLabel -> {"Edge Safety Factor ((Subscript[q, a]))","ω, γ"}, PlotLabel -> "External Kink for Wesson Profile, (Subscript[q, 0]) = (1.5, 1.4, 1.3)", PlotRange -> {0,20000}]

Plot[{Re[ωMHD[qa,1.5,1.8,6,2]],Im[ωMHD[qa,1.5,1.8,6,2]]}, {qa, 2.1, 3.2}];
Plot[{Re[ωMHD[qa,1.4,1.8,6,2]],Im[ωMHD[qa,1.4,1.8,6,2]]}, {qa, 2.1, 3.2}];
Plot[{Re[ωMHD[qa,1.3,1.8,6,2]],Im[ωMHD[qa,1.3,1.8,6,2]]}, {qa, 2.1, 3.2}];
Show[%%%,%%,%, Frame -> True, FrameLabel -> {"Edge Safety Factor ((Subscript[q, a]))","ω, γ"}, PlotLabel -> "External Kink for Wesson Profile, (Subscript[q, 0]) = (1.5, 1.4, 1.3)", PlotRange -> {0,20000}]

Internal Tearing Modes

Step 1: Define the equilibrium and the resonant surface

exampleRules = {qa -> 2.9, q0 -> 1.1, m -> 3, n -> 2}

q[x]

q[x] /. νj -> qa/q0 - 1 /. exampleRules

rs = x /. FindRoot[Evaluate[q[x]/. νj -> qa/q0 - 1 /. exampleRules] == 3/2, {x, 0.5}]

Plot[{q[x]/. νj -> qa/q0 - 1 /. exampleRules, 3/2}, {x, 0, 1}]

Step 2 : Find the inner solutions

eqsIdeal

origin

solInner = NDSolve[Evaluate[eqsIdeal ~Join~ origin /. νj -> qa/q0 - 1 /. exampleRules], {ψ[x], ψp[x]}, {x, δ0, rs-δ0}] // Flatten

ψrs = ψ[x] /. solInner /. x -> rs - δ0
ψp[x] /. solInner /. x -> rs - δ0

Plot[ψ[x] /. solInner , {x, δ0, rs-δ0}]
Plot[ψp[x] /. solInner , {x, δ0, rs-δ0}]

Step 2 : Find the Outer solution

bWall = 1.8;

outerBoundaries = {ψ[rs+δ0]==ψrs,ψp[1-δ0]==ψ[1-δ0] ψpψBoundary[1,bWall,3]}

solOuter = NDSolve[Evaluate[eqsIdeal ~Join~ outerBoundaries /. νj -> qa/q0 - 1 /. exampleRules], {ψ[x], ψp[x]}, {x,rs+δ0,1-δ0}] // Flatten

ψ[x] /. solOuter /. x -> rs + δ0
ψp[x] /. solOuter /. x -> rs + δ0

Plot[ψ[x] /. solOuter , {x, rs+δ0,1-δ0}]
Plot[ψp[x] /. solOuter , {x, rs+δ0,1-δ0}]

Step 3 : Show Full Solution

Evaluate[ψVac[xx] /. ψa -> ψ[x]/ψrs /. solOuter /. x -> 1-δ0 /. a -> 1 /. b -> bWall]

Plot[ψVac[xx] /. ψa -> ψ[x]/ψrs /. solOuter /. x -> 1-δ0 /. a -> 1 /. b -> bWall /. m -> 3, {xx,1, bWall}]

Plot[ψ[x]/ψrs /. solOuter , {x, rs+δ0,1-δ0}];
Plot[ψ[x]/ψrs /. solInner , {x, δ0, rs-δ0}];
Plot[ψVac[xx] /. ψa -> ψ[x]/ψrs /. solOuter /. x -> 1-δ0 /. a -> 1 /. b -> bWall/. m-> 3, {xx,1, bWall}];
Show[%%%,%%,%, PlotRange -> {{0,bWall},All}, Frame -> True, AxesOrigin -> {rs,1}]

Step 5 : Compute Δ'

ΔPrime = ((ψp[x] /. solOuter /. x -> rs + δ0) - (ψp[x] /. solInner /. x -> rs - δ0))/ψrs

Growth rate (with W/a ~ 0.05)

ΔPrime (η/μ0)/((a0^2) 0.05) /. hbtepRule
1/%

Plot[((ψp[x] /. solOuter /. x -> rs + w) - (ψp[x] /. solInner /. x -> rs - w))/ψrs, {w, δ0, 10 δ0}, PlotLabel -> "Δ' vs. W", AxesLabel -> {"W","Δ'"}]

More Tearing Modes

Step 1: Define the equilibrium and the resonant surface

mExample = 2;

exampleRules = {qa -> 4.2, q0 -> 1.1, m -> mExample, n -> 1}

q[x]

q[x] /. νj -> qa/q0 - 1 /. exampleRules

rs = x /. FindRoot[Evaluate[q[x]/. νj -> qa/q0 - 1 /. exampleRules] == Evaluate[m/n /. exampleRules], {x, 0.5}]

Plot[{q[x]/. νj -> qa/q0 - 1 /. exampleRules,m/n/. exampleRules}, {x, 0, 1}]

Step 2 : Find the inner solutions

eqsIdeal

origin

solInner = NDSolve[Evaluate[eqsIdeal ~Join~ origin /. νj -> qa/q0 - 1 /. exampleRules], {ψ[x], ψp[x]}, {x, δ0, rs-δ0}] // Flatten

ψrs = ψ[x] /. solInner /. x -> rs - δ0
ψp[x] /. solInner /. x -> rs - δ0

Plot[ψ[x] /. solInner , {x, δ0, rs-δ0}]
Plot[ψp[x] /. solInner , {x, δ0, rs-δ0}]

Step 2 : Find the Outer solution

bWall = 1.8;

outerBoundaries = {ψ[rs+δ0]==ψrs,ψp[1-δ0]==ψ[1-δ0] ψpψBoundary[1,bWall,mExample]}

solOuter = NDSolve[Evaluate[eqsIdeal ~Join~ outerBoundaries /. νj -> qa/q0 - 1 /. exampleRules], {ψ[x], ψp[x]}, {x,rs+δ0,1-δ0}] // Flatten

ψ[x] /. solOuter /. x -> rs + δ0
ψp[x] /. solOuter /. x -> rs + δ0

Plot[ψ[x] /. solOuter , {x, rs+δ0,1-δ0}]
Plot[ψp[x] /. solOuter , {x, rs+δ0,1-δ0}]

Step 3 : Show Full Solution

Evaluate[ψVac[xx] /. ψa -> ψ[x]/ψrs /. solOuter /. x -> 1-δ0 /. a -> 1 /. b -> bWall]

Plot[ψVac[xx] /. ψa -> ψ[x]/ψrs /. solOuter /. x -> 1-δ0 /. a -> 1 /. b -> bWall /.exampleRules, {xx,1, bWall}]

Plot[ψ[x]/ψrs /. solOuter , {x, rs+δ0,1-δ0}];
Plot[ψ[x]/ψrs /. solInner , {x, δ0, rs-δ0}];
Plot[ψVac[xx] /. ψa -> ψ[x]/ψrs /. solOuter /. x -> 1-δ0 /. a -> 1 /. b -> bWall/. exampleRules, {xx,1, bWall}];
Show[%%%,%%,%, PlotRange -> {{0,bWall},All}, Frame -> True, AxesOrigin -> {rs,1}]

Step 5 : Compute Δ'

ΔPrime = ((ψp[x] /. solOuter /. x -> rs + δ0) - (ψp[x] /. solInner /. x -> rs - δ0))/ψrs

Step 5 : Compute Δ'

ΔPrime = ((ψp[x] /. solOuter /. x -> rs + δ0) - (ψp[x] /. solInner /. x -> rs - δ0))/ψrs

Growth rate (with W/a ~ 0.05)

ΔPrime (η/μ0)/((a0^2) 0.05) /. hbtepRule
1/%

Plot[((ψp[x] /. solOuter /. x -> rs + w) - (ψp[x] /. solInner /. x -> rs - w))/ψrs, {w, δ0, 10 δ0}, PlotLabel -> "Δ' vs. W", AxesLabel -> {"W","Δ'"}]

Plot[((ψp[x] /. solOuter /. x -> rs + w) - (ψp[x] /. solInner /. x -> rs - w))/ψrs, {w, δ0, 10 δ0}, PlotLabel -> "Δ' vs. W", AxesLabel -> {"W","Δ'"}]

Summary

Cylindrical reduced MHD gives a relatively straight-forward method to compute kink and tearing instabilities in current-carrying tokamaks.